Temperature-dependent quantum pair potentials and their application to dense 

partially ionized hydrogen plasmas 
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Extending our previous work Q we present a detailed discussion of accuracy and practical ap- 
plications of finite-temperature pseudopotentials for two-component Coulomb systems. Different 
pseudopotentials are discussed: i) the diagonal Kelbg potential, ii) the off-diagonal Kelbg potential 
iii) the improved diagonal Kelbg potential, iv) an effective potential obtained with the Feynman- 
Kleinert variational principle v) the "exact" quantum pair potential derived from the two-particle 
density matrix. For the improved diagonal Kelbg potential a simple temperature dependent fit is 
derived which accurately reproduces the "exact" pair potential in the whole temperature range. 
The derived pseudopotentials are then used in path integral Monte Carlo (PIMC) and molecular 
dynamics (MD) simulations to obtain thermodynamical properties of strongly coupled hydrogen. It 
is demonstrated that classical MD simulations with spin-dependent interaction potentials for the 
electrons allow for an accurate description of the internal energy of hydrogen in the difficult regime 
of partial ionization down to the temperatures of about 60 000 K. Finally, we point out an interesting 
relation between the quantum potentials and effective potentials used in density functional theory. 



I. INTRODUCTION 

In recent years there is growing interest in the prop- 
erties of dense quantum plasmas, particulary in as- 
trophysics, laser plasmas, and condensed matter, see 
Refs. [SHSIEIEQl for an overview. In particular, the 
thermodynamic properties of hot dense plasmas are es- 
sential for the description of plasmas generated by strong 
lasers 0. Further, among the phenomena of current in- 
terest are the high-pressure compressibility of deuterium 
0, metallization of hydrogen and the h ypo thetical 
plasma phase transition, e.g. Ilfl[ Ei E El 13 E E3| 
which occur in situations where both interaction and 
quantum effects are relevant. 

While the case of strong degeneracy and the weak 
coupling limit have been extensively studied theoret- 
ically, e.g. within the random phase approximation, 
plasma properties at intermediate coupling and degen- 
eracy (when r - the ratio of the potential energy to the 
mean kinetic energy exceeds unity), are a hot topic of 
the present research activity. For an overview of present 
day analytical methods, see e.g. Refs. ^, ^, An- 
alytical methods typically use a chemical picture where 
electrons, ions and bound states (atoms, molecules etc.) 
are treated as independent species, and the chemical com- 
position (degree of ionization) is computed from a mass 
action law (non-ideal Saha equation). However, these 
methods are based on perturbation expansions in the 
coupling strength and are thus limited to regions of small 
coupling parameters, F < 1 or Ts < 1 (rs is the quantum 
coupling parameter, — f/as)- Furthermore, the mass 
action law becomes increasingly inaccurate in the region 
where the electrons are degenerate (uncertainty in the 
mass action constants). Also, during rapid pressure ion- 
ization around the Mott density the distinction between 



free and bound particles is an open problem. 

On the other hand, in the last decade, static proper- 
ties (e.g. equation of state) of dense hydrogen in ther- 
mal equilibrium have been successfully investigated with 
"exact" quantum-statistical methods, such as path inte- 
gral Monte Carlo (PIMC) [H El ElliOl- This first 
principle numerical technique is well suited for an ac- 
curate treatment of many-particle correlation effects in 
quantum systems, but unfortunately does not give dy- 
namical characteristics of the plasma (with an excep- 
tion of those obtained within the linear response the- 
ory). The alternative numerical approach for dense par- 
tially ionized plasmas (which does not have the above 
shortcoming) is a group of methods based upon ab ini- 
tio quasi-classical molecular dynamic simulations (MD), 
e.g. Refs. when a real quantum system is pro- 

jected onto a classical one where most of quantum ef- 
fects are included in some effective interparticle interac- 
tion potentials, such as the ones proposed by Kelbg [2fll |. 
Deutsch (2^, Klakow, Toepffer and Reinhard |2ll | and 
many others, e.g. UlElli^- These potentials can be 
derived from the two-particle Slater sum using Morita's 
method. 

However, no rigorous comparison of the accuracy of 
these potentials has been done yet, which is one of the 
aims of this paper. Different quantum potentials are 
compared with an "exact" pair potential obtained from 
the two-particle density matrix. Furthermore, we intro- 
duce pair potentials including particle statistics, e.g. de- 
scribing interaction between electrons in the singlet and 
triplet states, and use them in our MD simulations of 
two-component hydrogen plasmas. 

This paper is organized as follows: In Section |n] we 
discuss different methods to obtain an effective quantum 
pair potential. In the weak coupling limit, this poten- 
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tial leads exactly to the off-diagonal Kelbg potential the 
properties of which are being discussed and compared to 
its commonly used diagonal approximation. We outline 
two methods for a direct solution of the off-diagonal two- 
body Bloch equation which are then used in Section Hill 
for numerical comparison with the Kelbg and improved 
Kelbg potentials for rigorous assessment of the accuracy 
of the latter. Further, in Section HVl we present an analy- 
sis of the accuracy of the diagonal and off-diagonal Kelbg 
potentials in the PIMC simulations. Section Fvl describes 
an application of the improved Kelbg potentials to clas- 
sical molecular dynamics simulations of dense hydrogen. 
Comparing the results to those of PIMC simulations al- 
lows us to conclude that use of the improved Kelbg po- 
tential allows to significantly extend the range of appli- 
cability of classical MD to the region of partial ionization 
and to temperatures as low as approximately one third of 
the binding energy. Section discusses another field of 
potential applicability of the quantum potentials - den- 
sity functional theory. Finally, Section IVIII concludes the 
paper. 



II. EFFECTIVE QUANTUM PAIR 
POTENTIALS 



In this section we discuss different possibilities to ob- 
tain effective quantum potentials describing interactions 
in the two-particle problem. 



A. Analytical solution of two-body Bloch equation. 
Off-diagonal and diagonal Kelbg potential 

The equilibrium pair density matrix at a given in- 
verse temperature /3 = l/ksT is the solution of the two- 
particle Bloch equation 



— p(r,,rj,r-,r^;/3) = -iJ p(ri, , r r^; /3), 



H = k, + k,+U{r,,rj,r',,r'). 



(1) 



Numerical methods to obtain the density matrix of the 
Eq. 10 will be considered in Sec. Ill CI Here, we concen- 
trate on the available analytical solutions in the limit of 
weak coupHng. If the interaction is weak, the Eq. Q 
can be solved by perturbation theory with the following 
representation for the two-particle density matrix 



exp 



X exp 



exp[-/3$y], (2) 



wherei, j are particle indices, Pij = p{ri,rj,r[,r'j;P), and 
^ij = $(ri, Tj, r^, r' ; /3) is the off-diagonal two-particle 



effective potential. In the following we will consider ap- 
plication of this result to Coulomb systems. As a result 
of first-order perturbation theory we get explicitly 



<i>0(r. 



da 



-erf 



dij[a)/\ 



dij [a) \ 2 - a) 



(3) 



Q!)r'j |, erf(a;) is the error 



where dij [a) = \avij + {1 
function, erf(a;) — ^ dte~*^ , and Xfj = with 
The diagonal element (r^. = r^) 



TO- 



m. 



of (EJ is the potential derived by Kelbg and his cowork- 
ersUfSi 



VTTXy [1 - erf(a;ij)]| (4) 



with Xij = \rij\/Xij. The Kelbg potential is finite at zero 
distance reflecting that it captures the basic quantum 
diffraction effects and the quantum nature of two-particle 
interaction at small distances which prevents any diver- 
gence. From Eq. Q it is also clear that quantum effects 
become dominant (and there the quantum potential de- 
viates from the classical Coulomb potential) at distances 
Tij < Xij given by the thermal DeBroglie wavelength. 
We will see below that, in interacting systems, this is 
only a rough approximation, and at strong coupling, the 
expression for the quantum particle "extension" deviates 
strongly from Xij and needs to be generalized. 

To obtain a simplifled expression for the rather com- 
plex quantum potential lO one can approximate the off- 
diagonal matrix elements by the diagonal ones. A first 
possibility is to approximate the integral over a by the 
length of the interval multiplied with the integrand in 
the center (Mittelwertsatz) which leads to the so-called 
KTR-potential due to Klakow, Toepffer and Reinhard 
which (in the diagonal approximation) is often used in 
quasi-classical MD simulations 



1^<lj . ( dij{ll2) \ 

d.,(i/2) A., ; 



(5) 



where dij {1/2) 



2 1'-y 



Alternatively, the inte- 



gral can be simplified by taking the off-diagonal Kelbg 
potential only at the center coordinate. 



<i>°(r,„r^ /?)«$o. 



,/3 



(6) 



Many authors use the end-point approximation Q for 
the effective potential $(ry,r-^, /?) in the pair density 
matrix ^ due to the fact that it is very convenient com- 
putationally. The pair potential for interparticle inter- 
action simply is replaced by an effective potential which 
has only a dependence on the radial variables |ry |, |r-j|. 
However, most of the accuracy is usually lost in this end- 
point approximation. 
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Since the Kelbg potential is obtained by first order per- 
turbation theory its appHcation is hmited to weak cou- 
phng, r < 1, where F is the ratio of mean potential 
to kinetic energy. In unbound and bound states of an 
electron-proton pair this results in the following condi- 
tions on temperature 



-IksT < 1 



r = Ry/fcBr<l ^ fcBr>Ry, (7) 

where Ry = Ha/2 — e^/2aB, and gb is the Bohr radius. 
For the last case the Kelbg potential (and any of the 
simplifying approximations) can be only valid for tem- 
peratures sufficiently above the atomic binding energy, 
i.e for the case of hydrogen, T > Ry/fcs « 158 000 K. 
We address this point in more detail in the section ITVl 
where "exact" binding energies and pair correlation func- 
tions for an electron-proton pair are compared with the 
results obtained with the potentials (PJ and Q. 

B. Improved diagonal Kelbg potential 

The fimitation of the Kelbg potential to describe quan- 
tum systems only when there are no bound states has 
lead several researchers [ll l3(il8ll | to introduce and inves- 
tigate a more generalized form of the quantum potential 
with an additional free parameter 7^ 

^{nj,P) = (8) 



qzqj 



■•3 + 



1 - erf 



7y 



A. 



This potential has an advantage of preserving the correct 
first derivative at r = of the original Kelbg potential, 
$(0, PYj. = —qiqj/Xfj, but at the same time allows to 
correct the wrong value of the height of the Kelbg poten- 
tial at r = 0, i.e. $(0, (3) = qiqj\/T^ / jij) to include 
bound states. Using the definition of the effective poten- 
tial as 



(9) 

where Sij = S{rij,f3) is the exact binary Slater sum of 
particles The fit parameter 7^ in Eq. is related 
to the Slater sum at zero interparticle distance according 
to 



7ij 



(10) 



X^j ln[Stj{rij = 0,(3)] ' 

It is important to note that 7^ depends both on temper- 
ature and the type of particles. For example, the binary 
Slater sum of two electrons at zero separation has the 
form (including the average "(. . .)" over possible values 
of the total spin 5 = 0,1) 



5£(ree =0,/3) =2V^feeJl(eee), 



Jl(e., 



: dx 



l-exp(-^) 



(11) 



where the interaction parameter = qiqj P/Xij. 

On the other hand, for an electron-proton pair the 
Slater sum can be written as 



Sepirep = 0, /?) = 4V^eep Jl(Cep) + T^Cep ^3 (Cep ) , 



(12) 



where the last term shows the contribution of the bound 
states. 

The original Kelbg potential was derived for very high 
temperatures without taking into account exchange be- 
tween particles. This work was followed by several 
studies where the pseudopotentials for identical parti- 
cles have been calculated numerically or analyt- 
ically "s?, 35, "s^ using expansions in a quantum param- 
eter, small particle separation, and temperature. In the 
present work, following these studies, we approach the 
problem of the pseudopotential with exchange by using 
the formalism of two-particle density matrices (DM) . The 
pair DM can be calculated numerically (see Sec. lIICII or 
expressed in analytical form (01 using the improved Kelbg 
potential (jHJ. 

In the case of a pair of electrons, they can be in a sin- 
glet or triplet state, and the spatial wave function is sym- 
metric or antisymmetric under the exchange of particle 
indices. Thus, one can define a binary effective electron- 
electron interaction for three different cases 



S(T) 



{vi , Tj , r, , Tj ; /3) ± (r, , , , ; ^) 



3 -^^5 + 1 

4 4 



pW(r„r,;/3)pW(r„r,) 

(13) 



where p^^^ and p'^' are the one- and two-particle density 

matrices, and Uf,, Ujl and C/,-] are the effective interac- 

f-j ''J 

tions in the singlet (S), triplet (T) state and the spin- 
averaged potential, respectively. 

If we now approximate the two-particle DM, p^^\ by 
Eq. ((2l and factorize it into the DM's of the center-of- 
mass and relative coordinates (the corresponding expres- 
sions are given in Sec. Ill ("I cf. Eq. Il27|l l. then we obtain 
for the pseudopotential between two electrons being in 
the singlet (triplet) state and for the spin-averaged po- 
tential, respectively, 

UO = _iin (^e-'^^-f--^--) - ^e-'-'/^" e-'3t^"(--'--))(15) 

In this expression, the function, J7ee(r,r'), is a pseu- 
dopential between distinguishable particles (i.e. calcu- 
lated without exchange effects). Thus, for it one can 
substitute the original Kelbg potential, Eq. 1^, the im- 
proved Kelbg potential, Eq. (|Hl, or any further improved 
approximation for the binary interaction. In the case of 
two electrons, if for J7ee(r,r) we use the improved Kelbg 
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potential Q, then the fit parameter 7ee must be ob- 
tained from Eq. ifTTHl where for the binary Slater sum 
one should take the one for two distinguishable particles 
with Coulomb repulsion 

^e"°'"'(^ee = 0, /?) = ^V^^eeMLe)- (16) 

As it follows from Eq. lfT5)l an exchange contribution 
(effect of particle statistics in the pair interaction) arises 
from the kinetic energy part of the density matrix and 
the non-diagonal potential, J7ee(r, — r), which in the first 
order of the perturbation theory can be calculated us- 
ing Eq. (PJ. A further simplification (which is crucial 
for application of the pseudopotentials in semiclassical 
MD simulations presented in Sec. can be achieved 
by approximating the off-diagonal potential by the di- 
agonal terms, f/ee(r, -r)« i [C/ee(r,r) -I- C/ee(-r, -r)] = 
?7ee(r,r), and the above expressions are reduced to 



3.0 




S(T) 



ee,0 



u. 



= C/ee(r,r)-iln{l±e-'-'/^==}, (17) 

(18) 



ee,o ^ C^ee(r,r) - -In <{ 1 



We can note that in the diagonal approximation for the 
potential, the exchange term corresponds to the case of 
the ideal Fermi gas (i.e. exchange without interaction), 
the exchange term arising from the interaction is missing. 

Taking in Eq. ifTTjl the limit r ^ we see that the po- 
tential of the triplet state shows a logarithmic divergency 



C/J,_o = C/ee(r,r)-2fcsTln 



O 



r 

A?, 



(19) 



whereas the singlet and the spin-averaged potential ac- 
quires an additional exchange contribution 



u: 



s,{) _ 



Uee{r,r)TkBT\n{2} 



(20) 



but the slope of these potentials at the origin are same 
as in the case without exchange. This means, in case of 
Coulomb interaction, the slope is defined by the slope of 
the original Kelbg potential 



C/ee(r,r)|r^O 



$°e(0)-^ 



(21) 



In our previous paper Q we reported on the tem- 
perature dependence of the fitting parameter jij for 
the electron-electron and electron-proton interactions. 
There, two types of calculations have been presented. 
First, the values of 7(/3) were obtained by a least-square 
fit of the improved diagonal Kelbg potential (IDKP), the 
Eq. l|Sl, to the "exact" pair potential U (see Eq. l(27jl l. and 
second, from Eq. lfT?Hl by evaluating the values of the bi- 
nary Slater sums. It has been found that both methods 
agree within the statistical uncertainty. 



Figure 1: Temperature dependence of the fit parameter for 
the binary interactions: electron-proton, 7ep(r), and electron- 
electron (no exchange), "fee{T). Symbols show the 7- values 
obtained with the least-square fit of the IDKP to the "exact" 
pair potential without exchange. Solid curves correspond to 
the Fade approximation 12212311 . 



Extending our earlier results, we now present a Fade 
approximation which contains an analytical temperature 
dependence of the parameters jij which will be useful for 
practical applications. 



lee{T) 



Xi + xf 



1 + aep xi + xi 



7ee(T ^ 0) + flee Xx ^ : 



(22) 

(23) 



where xi — -^/STrfcsT/Ha (with the Hartree energy. 
Ha = 2Ry = 315 775 K), a^p = 1.090(14), Oee = 0.18(1). 
The limit value, 7ee(T 0), has been obtained from 
Eq. lfT?Hl by evaluating the zero temperature limit of the 
binary Slater sum 



2 ^3 1 
^"""T^"" ln{8xy^} 



3x 



~2 ■ 



(24) 



with X = (|7r^ee|/2)^. The excellent accuracy of the Fade 
approximation is demonstrated in Fig. 

In Fig. Q we present the temperature dependence of 
7 obtained from the least-square fit (full and open sym- 
bols) to the "exact" pair potential of distinguishable par- 
ticles (no exchange) and Fade approximation II22II (solid 
curves). The most important result is that the corrected 
Kelbg potential is now not limited to weak coupling as 
is the original Kelbg potential. For the case jij — 1, 
Eq. (HJ coincides with Eq. Q. One clearly sees the de- 
viation of 7ij from unity for T < 10^ K, which shows 
that the quantum extension of particles is becoming in- 
fiuenced by interaction effects and is now of the order of 
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Aij = ^ijlij^ instead of the original thermal DeBroglie 
wavelength Ay . Thus, with the Fade formulas H22I23II . 
we have obtained an analytical fit for the quantum ex- 
tension of the scattering particles. 

The Fade approximations lf22ll , ll23ll have been success- 
fully used in quasi-classical molecular dynamics simula- 
tions of two-component hydrogen plasmas. As we show 
in Sec. they allow to obtain accurate results for par- 
tially hydrogen (or other quantum systems of oppositely 
charged particles with bound states). 

Finally, we want to note that, in the MD simulations of 
quantum plasmas, it can be advantageous to have spin- 
dependent potentials for the electron subsystem defined 
by Eqs. ^ ov The spin-resolved approximation 

allows for a refined modelling, for example, it allows for 
the description of molecule formation, spin density waves, 
spin fiip processes in the presence of a magnetic field and 
so on. 



C. Effective potentials from numerical solution of 
two-body Bloch equation 

In the present section, we briefiy describe the numerical 
methods which have been used to solve the two-particle 
problem in order to obtain the "exact" quantum pair po- 
tentials. These results have been used to obtain the ana- 
lytical fit in the improved Kelbg potential. Furthermore, 
they will be used to test the accuracy of various analyt- 
ical approximations for the quantum pair potentials in 
Sec, mil below. 

Let us factorize full two-particle density matrix into 
a center-of-mass term and a density matrix of relative 
coordinates 

p(r„r,,r^,r;;/3) =pcm(R,R';/3)p(r,r';/3), (25) 

where R = {rmvi + rajYj)/{mi + raj), and r = — r^, 
and analogously for R', r'. When for the relative DM in 
the analogy with the Eq. Q we can define the effective 
pair potential as 

p(r, r'; P) = pfc,„(r, r'; /3) e''^^^'^-'-'^'^), (26) 

which results in the following expression 

C/(r, r'; /3) = - ^ In [p(r, r'; /3)/pfc„ (r, r'; /?)] . (27) 

where p(r, r'; I3)kin is the kinetic energy DM. 

One of the possibilities to get the relative density ma- 
trix p(r, r'; (3) is to directly solve the corresponding one- 
particle Schrodinger equation and calculate the DM as 
a contribution from bound and continuum states. This 
procedure is advantageous when the Schrodinger equa- 
tion can be solved analytically and we know analyti- 
cal expressions for contributions of scattering and bound 
states, as for example for Coulomb potential, e.g. [3^ . 
But if it is not the separate calculation of each 



matrix element for each new values of end-points r and 
r' will be not efficient and time-demanding procedure. In 
principle, such calculations can be done in advance with 
C/(r,r';/3) be stored in the tables of the potential, but 
one still needs to solve the Schrodinger equation many 
times for each value of quantum numbers and also for 
wavefunctions of continuum states. 

It is possible to approach to this problem from the 
other side and calculate the DM directly without solving 
the Schrodinger equation. In this work we apply two 
efficient methods - the matrix squaring technique 
l^a and the Feynman-Kleinert variational approach [3^ 
l3ijl |. In the Sec. IIIII we will compare an accuracy of the 
pseudopotnetials obtained with these methods. 



1. Matrix squaring technique 

The exact off-diagonal pair density matrix can be cal- 
culated efficiently by this method introduced by Storer 
and Klemm j^^l- For the case of spherical symmetry of 
the interaction potential, the relative pair density matrix 
in the Eq. Ij25|l is expanded in terms of partial waves. 
This expansion reads, for the two- and three-dimensional 
cases. 



+00 



p2^(r,r';/3) = ^ pi{r,r' ■ P)e'''' , (28) 

ZTT V r r — 



/ — — oo 

+ 00 



p3^(r,r';/3) = ^-^y^ (21 + I) pi{Ty ■ P) Pi{cosQ), 
Anr r ^-^ 

1=0 

where 6 is the angle between r and r'. Each partial- wave 
component satisfies the ID Bloch equation for a single 
particle in an external potential given by the interaction 
potential and also a convolution equation, 

oo 

pi{r,r';T)^ f dr" pi{r,r" ;t/2) piir" ,r' ■,t/2). (29) 



This is the basic equation of the matrix- squaring method 
which allows to calculate the function pi at a given tem- 
perature l/r from the same function at a two times 
higher temperature. Squaring the density matrix k times 
results in a lowering of the temperature by a factor of 
2'^. Each squaring involves only a one-dimensional in- 
tegration which, due to the Gaussian-Hke nature of the 
integrand in Eq. l|29ll . can be performed quite accurately 
and efficiently by standard numerical procedures. To 
start the matrix-squaring iterations, Eq. ll29ll . one needs 
a known accurate high-temperature form for the density 
matrix. A convenient choice is the semiclassical approx- 
imation 



pi{r,r']T) = p^{r,r';T)exp - 



\r — r'\ 



V{x)dx 



(30) 
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where p'^{r,r';T) is the partial- wave component of the 
free-particle density matrix. 

Once the pair density matrix pi{r,r';T) is computed 
for the desired value of r, it is substituted into Eqs. lf28l 
I29II . and a summation over partial waves readily yields 
the full relative density matrix. 



2. Variational perturbation approach 

As a second method for solving the off-diagonal Bloch 
equation we used a variational perturbation expansion de- 
veloped by Feynman and Kleinert In this procedure 
the initial density matrix is presented in the form of a 
trial path integral which consists of a suitable superposi- 
tion of local harmonic oscillator path integrals centered 
at arbitrary average positions Xm , each with its own fre- 
quency squared ri^(xm). One starts from decomposing 
the action in the density matrix as 

p(r,r';/?) = J I?xe-^W/^ (31) 

(r,0)^(r',ft/3) 

^[x] = Ao,x„[x] -I- A„t[x], (32) 

with Afi^ x„ [x] being the action of a trial harmonic os- 
cillator with the potential minimum located at Xm, and 
V being the functional integral over all trajectories. The 
interaction part 



Aj«t[x] = / ' 
"'0 



drj 



F[x(r?)]- 1^172 [x(7;)-xj2 



(33) 

is defined as the difference between the original potential 
y(x) and the displaced harmonic oscillator. The fi^ term 
in Eq. (l33ll compensates the contribution of Aq^ x„ [x] in 
Eq. l(32|l . Now one can calculate the density matrix l(3T|l 
by treating the interaction (j33ll as a perturbation, leading 
to a moment expansion 

p(r,r';/3) = (r, r'; /3) (^1 - i(A,„,[x])^r''" + 



d/2 



(34) 



with the definition 



0,x„ 
N 



d J sinh h|3^l 
X [(r^ +y''^) coshhpn - 2rr'] - 

-^Ey?(^»*M)r'!r'^ (35) 
n— 1 



where d is the space dimensionality and A'' the order of 
the approximation. The function pg ' (r, r') is the trial 



harmonic oscillator density matrix, r = (r — Xm), ?' = 
(r' — Xm), and the expectation value of the interaction 
action on the r.h.s. of Eq. (l35ll is given by 



j V~.\{[jdn. 

'.V^^t [x(rO +x„] e{-*^n.==™[x+x„,]}|^ (3g^ 



The function W] 



n, x„ 

N 



can be identified as an effective 
quantum potential which is to be optimized with respect 
to the variational parameters {J7^(r, r'; /3), Xm(r, r'; /?)}. 
Note that, in the high temperature limit, this effective 
potential goes over to the original potential ^(r). The 
optimal parameter values are determined from the ex- 
tremity conditions 



5<-""'(r,r^/?) 



0, 



9<-"'"(r,r';/3) 

9Xm 



0. (37) 



The perturbation series Ij35ll is rapidly converging, in 
most cases already in the first-order approximation 
W^' ^" for the effective potential, and gives a reasonable 
estimate of the desired quantities. 



III. COMPARISON OF THE PAIR 
POTENTIALS AND THEIR TEMPERATURE 
DEPENDENCE 

We will now compare accuracy of the pair potentials 
discussed above (or two-particle density matrices corre- 
sponding to these potentials), their temperature depen- 
dence and range of applicability. 



A. Full density matrix of electron-proton pair 

In Fig.|2lwe show the angular dependence of the full off- 
diagonal two-particle density matrix calculated with the 
off-diagonal Kelbg potential - ODKP ^ and its diagonal 
approximation - DKP I^J • The density matrix is shown 
at several temperature values (T = 1 000 000, 250 000 
and 62 500 K) and several angular distances (0 — 
0, 7r/2, tt) between the vectors r = rtj, r' = r'^j (in each 
of the figures, the top curves correspond to the case of 
parallel vectors, (j) — 0, the lowest curves to antipar- 
allel vectors, cj) — n). Also, for reference, we give the 
off-diagonal density matrix obtained from the 'exact' so- 



lution of the Bloch equation, cf. Sec. Ill C II At high 
temperatures, T > 250 000 K, the Kelbg density matrix 
does not exhibit large deviations from the exact result. 
At T = 1 000 OOOK, the ODKP density matrix practi- 
cally coincides with the exact solution, whereas the DKP 
approximation shows small deviations. In these cases the 
perturbation expansion applies, F ~ 0.15, see left column 
of Fig. 12 With decreasing temperature, the deviations 
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Figure 2: The "exact" off-diagonal density matrix p{r,r'\<f>) 
for an electron-proton pair vs. the density matrix calcu- 
lated with the diagonal (DKP) and off-diagonal (ODKP) 
Kelbg potentials. In all figures, results for three angular 
values are given = (upper curves), = 7r/2 (middle) 
and <f> — -K (lower curves). The proton is located at the 
origin, and the vector r' (initial electron position) is fixed, 
|r'| = 0.25; 1.0; 2.0. The vector r (final electron position) is 
varied, 4> is the angle between the vectors r and r'. 



from the exact results grow, see middle column. To bet- 
ter understand the details of the deviations, we magni- 
fied them by including also results for T — &2 500 K, 
which is far beyond the scope of the perturbation theory, 
T « 0.4 Ry/fcs, i.e. F « 2.5. Here we observe that, at 
the origin, the density matrix of the Kelbg potential is 
3 times less than the exact one. The largest errors were 
found for the DKP, in particular, in the case when the 
vectors r,r' have the opposite direction ((/) = tt). 

The behavior of the full density matrix can be under- 
stood from the following considerations. The density ma- 
trix results from contributions of kinetic and potential en- 
ergy operators, cf. Eq. QJ. At small distances (r' = 0.25) 
the Coulomb attraction between an electron and a pro- 
ton dominates and, therefore, the density matrix shows 
an exponential decay. At the largest distance (r' = 2.0) 
kinetic and potential energy are of the same order and 
a Gaussian-like free particle density matrix emerges, as 
can be clearly seen in the bottom left part of Fig. [2 

From this first comparison we can conclude that both, 
the DKP and the ODKP, show satisfactory agreement 
with the exact result in the cases where perturbation 
theory applies, T > 2Ry. At lower temperatures there 
is only qualitative agreement. The strongest deviations 
arise for small interparticle distances {r, r'}, and this, as 



will be shown below, results from the incorrect height of 
the Kelbg potential at zero distance r = 0. 



B. Effective interaction of electron-proton and 
electron-electron pairs 

In Fig. |2ta,b) we show and compare the accuracy of 
several effective electron-proton potentials and their tem- 
perature derivatives obtained by various methods. As an 
"exact" reference potential to which the accuracy of other 
potentials is compared we use Upair obtained from the 
electron-proton pair density matrix calculated with the 
matrix squaring technique. 

First, we note from Fig. EJa) that, at given tempera- 
tures T < 2Ry, the original Kelbg potential shows the 
largest deviations from the "exact" result, Upair- While 
the spatial derivative of the DKP coincides with that of 
Upair, a systematic offset of the DKP compared to Upair 
is observed at the origin r = 0, which increases when the 
temperature is lowered. The agreement is satisfactory 
only for the curve corresponding to T = 320 000 K. The 
accuracy of the Kelbg potential becomes worse for quan- 
tities involving its temperature derivative. For exam- 
ple, for the total energy one has to compute the thermo- 
dynamic average of the function d{P U{f3)/d(3)\u^^o^ . 
This function is shown in Fig.E^b). If multiplied by the 
Boltzmann factor e~^'^^^\ this function is a good esti- 
mate for the binding energy (Eb) of an electron-proton 
pair. In the case of a bound state the main contribu- 
tion to the energy comes from the region of small inter- 
particle distance, r < 3aB. Therefore, the behavior of 
d{(3 U{l3)/d(3) near the origin determines the accuracy of 
the calculations of the energy and other thermodynamic 
quantities. As we can see from the curve for 5 000if in 
Fig- 13b), the depth of the DKP is much less than that 
of Upair and, therefore, it gives a too low binding energy 
of Eb « 0.16 Ha, i.e. a factor of three too low compared 
with the true ground state energy E^ = 0.5 Ha. 

As it was already discussed in Sec. IIIBI the accuracy 
of the DKP can be improved with the additional fit pa- 
rameter 7y . In Fig. |Hta,b) this potential is denoted as 
Kelbg- One Can see that at all considered temperatures 
^ Kelbg practically coincides with Upair - Even in the case 
of the strong coupling (T = 5 000 K) the agreement is 
very good. 

The next potential showh in this figure is the varia- 
tional potential, W'^'^'"-, introduced in the Sec. Ill C 2l 
This potential is more accurate than the DKP and qual- 
itatively reproduces the "exact" effective pair potential 
Upair for temperatures T = 125 000 K, 320 000 K and its 
derivative (see Fig.EJb)). The key point is that the vari- 
ational perturbation theory [38j replaces the perturbation 
expansion in F (which does not converge for F > 1) into 
another expansion, Eq. Il34|l . which does not have this re- 
striction. The results of this approach can be improved 
by taking into account higher order terms in Eq. l(3ijl (the 
results shown in the figure include only the first term. 
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n = 1). The convergence of Eq. lI'MII extends even to 
very strong coupling and has been successfully appHed in 
field theory 0. 




Figure 3: (a): Diagonal effective electron-proton potential 
(in units of Ha) for several cases: the DKP $°(r;/3) the 
improved DKP $(r; (i) JSJ, variational potential W^^'^™ 
pair potential Up \27\ corresponding to the "exact" density 
matrix. Each potential is given at at three temperature val- 
ues: 5 000,40 000,125 000 and 320000 K. (b): the function 
U{P) + pdU{P) / dp for the same approximations and temper- 
atures. 

We mention that comparison with other effective po- 
tentials has been performed in our paper In partic- 
ular, the "exact" pair potential was compared with the 
results of Barker (the calculations of the pair po- 
tential by the direct eigenfunction summation) and the 
Deutsch potential. A good agreement has been found 
with the data of , while the deviations of the Deutsch 
potential turned out to be sHghtly larger than that of the 
Kelbg potential. The reason is that the Deutsch potential 
has an incorrect spatial derivative, U'^, for r < Sos. 

Next, in Fig. 31 we compare pair distribution func- 
tions (PDF) of two electrons in a singlet and triplet states 
for different temperatures, obtained from the expression 
with the effective potential 

5(r)cxe-^^''"'«. (38) 

Due to the Pauli principle in Fig. the PDF is goes 
to zero as Tee 0. On the other hand, for electrons in 
a singlet state (Fig. b) this happens only if the tempera- 
ture is decreased down to 31 250 K, when the potential 
energy dominates over the kinetic energy. The three lines 
in Fig.^^a) show the cases when as an effective potential 
in the Eqs. (I15I17|I we substitute the "exact" pair poten- 
tial and the Kelbg potential. In the last case, Eq. itTTIl . 
the exchange contribution from the potential function is 
neglected. This, as shows Fig. 0ta), becomes important 
only for the temperature 31 250 K and below. But the 
overall accuracy of the Kelbg potential for a description 
of two particles with the same charge (even without im- 
proving its value at the origin with 7) is significantly 
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Figure 4: Pair distribution function for a pair of electrons, a: 
In a triplet state (parallel spins, S = 1); b: In a singlet state 
(antiparallel spins, S — 0). Cases compared are: Eq. H14|l. us- 
ing for Uee the "exact" pair potential and the Kelbg potential, 
'I'ee.o, respectively, as well as the simpler expression, Eq. I|17|l 
with the Kelbg potential. The pair distribution functions are 
shown for four temperatures, T = 31 250, 125 000, 250 000 
and 1000 000 K. 



better (compared to the results with the "exact" pair po- 
tential) than for particles with opposite charge, cf. Fig.El 
This is due to absence, in this case, of contributions from 
bound states. 

The pair correlation functions at T = 125 000 K can be 
compared with those for a hydrogen plasma obtained by 
molecular dynamics simulations, see Fig.|2l One can see 
that, at least at this temperature, the electron-electron 
PDF's of that many-particle simulation look quite similar 
to the two-particle case shown in Fig. 31 Some differences 
are noticeable for the value of gee at r = as the density 
is increased from rs = 6 to Ts = 4, see Sec. 

In the next section we discuss application of quan- 
tum pair potentials in thermodynamical calculations us- 
ing Feynman trajectories in imaginary time (PIMC). 



IV. QUANTUM PAIR POTENTIALS IN PATH 
INTEGRAL MONTE CARLO 

It is well known (see for example discussion in Chapter 
12 of jl^) that the singularity of the attractive Coulomb 
potential causes difficulties in the euclidian path inte- 
gral. If based on Feynman's original path integral repre- 
sentation, a path consists of a finite number of straight 
pieces, each with a classical euclidian action, containing 
the singular Coulomb potential. However, in this case, 
the energy of the path can be lowered indefinitely by an 
almost stretched configuration which corresponds to a 
slowly moving particle sliding down to the — e^/r abyss. 
This phenomenon is called path collapse. 

One possibility to prevent this effect is to use a mod- 
ified "regularized" Coulomb potential which has a cutoff 
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at r = 0. This procedure, however, is quite arbitrary, 
and the results are sensitive to the used cutoff parame- 
ters. Of course, in nature, these difficulties are prevented 
by quantum fluctuations which equip the path with a 
conflgurational entropy. The latter must be sufficiently 
singular to produce a regular free energy bounded from 
below. The inclusion of quantum fluctuations in the Eu- 
clidean action of the Feynman path pieces smoothes the 
singular Coulomb potential, producing an effective po- 
tential that is flnite at the origin, and the path collapse is 
avoided. This again shows the importance of effective po- 
tentials, speciflcally, in "quasiclassical" simulations (clas- 
sical Monte Carlo and molecular dynamics methods). Of 
great importance are potentials which have a closed an- 
alytical form. In this case for many thermodynamical 
quantities it is possible to obtain analytical solutions. 

For simulations of correlated quantum many-body sys- 
tems which are based on first principles, the initial many- 
body hamiltonian with the true singular Coulomb energy 
operator is to be considered and solved to flnd some 
effective many-body interaction potential. For this it 
is important that in the high-temperature limit the N- 
particle density matrix can be expanded in terms of 2- 
particle, 3-particle etc. contributions. If the temperature 
is sufficiently high then all contributions except the first 
one, taking into account two-particle correlations, can be 
omitted. As a result, the following approximation for the 
N-particle density matrix holds 



N 



P(R,R';t) « l[p^'Hr^,r',■,r) 



pP](rj,rfe,r^.,r;^;T) 



)pW(rfe,r;;r) 



(39) 



0(p[3l), 



where R = {ri,...,r7v} specifies coordinates of all TV 
particles, p'^' (p^^') is the single (two) particle density 
matrix. The above pair approximation is usually used 
in PIMC simulations The N-particle density matrix 
p{(3) contains a complete information about the system 
with the observables given by 



(0) = 



Tr 



Trim] 



dR (R|C'/5(/3)|R) 
dR (R|p(/3)|R) 



(40) 



Due to the exponential form, the N-particle density op- 
erator p(/3) = e~^^ can be factorized (in analogy to the 
matrix squaring method above) as /5(/3) = [/5(r)]*^ with 
M = i3/t. Consequently, the N-particle density oper- 
ator p{(3) is expressed in terms of density operators at 
an M times higher temperature l/r = M ■ ksT . If M 
is chosen sufficiently large then one can apply the pair 
approximation Il39|l . Thus, accurate results for the quan- 
tum pair potentials and, consequently, the pair density 
matrix, will allow to compute the density matrix of the 
whole iV— particle system. Here we are not interested in 



the investigation of the accuracy of approximation Ij39|l 
but concentrate on the two-body problem where Eq. l(39jl 
is exact. 

It is clear that the observables Ij40ll computed with the 
approximate pair-density matrix pl^l contain an error of 
the order 0(1/M^). Below we will investigate the con- 
vergence, as a function of A/, of main thermodynamic 
properties (total energy and e-p pair distribution) for an 
electron-proton pair using for the pair density matrix p'^' 
results computed with the off-diagonal and the diagonal 
Kelbg potential. 



A. Comparison of diagonal and off-diagonal Kelbg 
potential on the example of Hydrogen atom 

We consider a hydrogen atom in a box with periodic 
boundary conditions (box size, L = 20 as) at several 
temperatures, T = 31 250 — 62 500 K, when the hydrogen 
atom can ionize into free particles, as well as for the case 
T < 10 000 K, when there is essentially only the contri- 
bution from the atomic ground state. First, in Fig. we 
show the e-p pair distribution functions (normalized to 
the volume dV = Airr'^dr). 
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Figure 5: Proton-electron pair correlation functions from 
PIMC simulations with the "exact" pair density matrix (di- 
amonds), the DKP (dots) and the ODKP (full line). Tem- 
perature values are as indicated in the figure parts: T — 
5 208, 10 416, 31250 and 62 500 K. For comparison, '3D GST' 
- denotes the pair correlation function corresponding to the 
ground state of a hydrogen atom. 

For temperatures T = 5 208 K and 10 416 K the hy- 
drogen atom does not decay into free particles during 
the duration of a typical simulation run (~ 10^ Monte 
Carlo steps). In the figures, the "exact" pair correlation 
function is compared with the one obtained with the off- 
diagonal and diagonal Kelbg potentials, respectively (the 
number of factorization factors for the density matrix was 
chosen to be M = 400). We found that the best accu- 
racy is achieved for the off-diagonal Kelbg potential and 
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M > 200, in this case the ODKP pair correlation func- 
tion is very close to the exact one. 

At elevated temperatures, T = 31 250 K and 62 500 K, 
ionization of the hydrogen atom occurs, but due to the 
periodic boundary conditions, the free particles cannot go 
to infinity but, when reaching the boundary, are returned 
back in the simulation box and have a finite probability 
for a formation of a bound state again. Thus, this sim- 
ulation captures the region of partial ionization. As the 
temperature is increased the ionization probability also 
increases, leading to a significant drop in the height of 
the proton-electron pair distribution function at the ori- 
gin compared to the ground state probability function 
^-§(0, (see Fig.El plots for T=31250 K and T=62 500 
K). 



62 500 K 1 41 6 K 




M M 



Figure 6: PIMC results for the internal energy of the proton- 
electron pair using "exact" pair density matrix, DKP, ODKP 
and Improved Kelbg potential vs. different number of factor- 
ization factors M. 



of the relative error, \og{SE/E), is shown as a func- 
tion of the inverse of the temperature used in the high- 
temperature factors, 1/r. In this figure we compare the 
behavior of the error for the same set of temperatures 
as in Fig. El In Fig. El we also add simulation results 
using improved diagonal Kelbg potential (solid Hue). Its 
accuracy is better then that of the ODKP at low tem- 
peratures (small values of M) but at high temperatures 
both are comparable. 

The main conclusion that can be drawn from the pre- 
sented PIMC results is that, at equal number of fac- 
torization factors M, simulations with the off-diagonal 
Kelbg potential are significantly more accurate in repro- 
ducing the "exact" thermodynamic results of a hydro- 
gen atom. Besides, the full off-diagonal density matrix 
contains valuable information about the spatial electron 
distribution around the proton, which is lost in the end- 
point approximation. Further, we expect that the best 
results will be obtained using an improved off-diagonal 
Kelbg potential, which has the correct zero-point value 
and contains the complete angular dependence of the pair 
density matrix which, however, is beyond the scope of the 
present paper. 



o DKP 
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In Fig. El we analyze the convergence of the internal en- 
ergy in PIMC simulations with varying number of high- 
temperature factors M. In particular, we compare inde- 
pendent simulations with the diagonal and off-diagonal 
Kelbg density matrices, respectively. The "exact" energy 
value for the considered temperatures is given by the soHd 
line and is obtained from PIMC simulations using the 
"exact" pair density matrix, cf. Sec. Ill ("Hi . The internal 
energy was obtained using the thermodynamic estimator, 
(E) = — ^ InZ, where Z is the partition function. Com- 
paring the diagonal and off-diagonal cases one can note 
that the ODKP density matrix shows much better and 
faster convergence to the exact energy value. A simple 
estimate shows that the relative error of the total energy, 
in the diagonal approximation, depends on factorization 
number M as SE/E « 7t^, t = j3/M. In contrast, us- 
ing the off-diagonal potential, the error converges much 
faster, SE/Ew-fT^. 

This fact is illustrated in Fig. where the logarithm 



Figure 7: Relative error in the internal energy of the proton- 
electron pair from PIMC simulations with diagonal, off- 
diagonal Kelbg potential vs. temperature argument in the 
two-particle density matrix 1/r. 



V. MOLECULAR DYNAMICS SIMULATIONS 

In this section, we apply the improved Kelbg poten- 
tials in classical molecular dynamics simulations of dense 
hydrogen. Classical MD simulations of dense plasmas 
have been performed by many authors before where the 
divergency at zero distance leading, in particular, to the 
classical collapse of an electron into a proton, is usually 
avoided by some cutoff or "regularization" of the Coulomb 
potential at small distances. Using in this paper, effec- 
tive quantum pair potentials obtained from the exact so- 
lutions of the Bloch equations, we expect to be on well 
founded grounds regarding the correct pair interactions 
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at short distances. This should not only prevent any 
collapse but also correctly reproduce the formation of 
hydrogen atoms and thus allow us to obtain essentially 
improved MD simulation results. On the other hand, this 
is not a trivial question, since these potentials are derived 
from pure equilibrium considerations, and there is no ad- 
hoc proof that they will necessarily be accurate for the 
description of dynamical behavior as well, in particular 
under strong non-equilibrium conditions. We, therefore, 
concentrate in the present analysis on correlated partially 
ionized hydrogen in thermodynamic equilibrium. The re- 
sults obtained below confirm that, indeed (at least in 
equilibrium) , the quantum pair potentials are well suited 
for use in the interparticle force terms in classical MD. 

Classical MD simulations incorporate all interparti- 
cle collision processes and are thus not restricted with 
respect to the coupHng parameter F in a classical sys- 
tem. With the use of effective quantum pair potentials, 
we expect, in addition, to capture dominant features 
of the quantum nature of microparticles, such as quan- 
tum diffraction and spin effects. Thus, these simulations 
could be called "semiclassical" MD. Having access not 
only to improved electron-ion potentials but also to spin- 
dependent electron-electron potentials, allows us to con- 
sider also spin effects by simulating electrons with differ- 
ent spin projections as two independent particle species. 
No spin-fiip processes, as they would be expected e.g. in 
strong magnetic fields are considered but our model 
should be capable to describe related phenomena as well. 
In this paper we focus on static properties, such as inter- 
nal energy, and radial distribution functions. Investiga- 
tion of dynamical properties an of spin density fiuctua- 
tion is the aim of a forthcoming paper. 

We consider a dense, degenerate hydrogen plasma at 
two densities corresponding to the Brueckner parame- 
ter rs = f/aB — 4 and Ts = 6 and temperatures 
T = 31 250, 50 000, 62 500, 125 000 and 166 670 K. These 
parameter values correspond, respectively, to r= 2.53, 
1.58, 1.26, 0.63 and 0.47, for r,= 4, and 1.68, 1.05, 0.84, 
0.42, 0.32, for r^= 6. 

The simulation box of our system, with the length 
L=[n/{Np + Nj + A^i)]^/^, contains Np= 200 protons, 
NI= 100 electrons with spin up and an equal number of 
electrons, N^= 100, with spin down. We keep the con- 
dition of the electro neutrality by taking Np= Nj. 
Details of the used numerical algorithm can be found in 
Ref. 

Since MD, in contrast to PIMC, involves only diag- 
onal interaction potentials, we choose the following ex- 
pressions: for the interaction between electrons and pro- 
tons, protons and protons and electrons with opposite 
spin, we use the improved Kelbg potential, Eq. {HI, with 
the fit parameters given by Eqs. ll22ll and ll23ll . respec- 
tively. The interaction between electrons with the same 
spin projection is described by the diagonal antisymmet- 
ric potential, Eq. ifTTjl . Further, to properly account for 
the long range-character of the potentials, we used the 
standard Ewald procedure as in Ref. [221 whereas, in 



contrast to the rather involved expressions given there 
for a one-component plasma, here we could restrict the 
potential energy sum only by the proper sum in real space 
(we do not reproduce these lengthy expressions here, but 
mention that the value of the parameter a defined in 
Ref. [131 was chosen to be a = 5.6/L) and taking 5 
vectors in every direction in the reciprocal space. This 
gives some computational-cost advantage in computation 
of the forces compared to Ref. (2^ . 

In Fig. IHl we plot the internal energy per atom as a 
function of temperature for two densities and compare it 
to the path integral Monte Carlo results of Militzer [iH . 
One can note, that for high temperatures the energies 
of MD and PIMC simulations coincide very well and He 
within the limits of the statistical errors. This is an im- 
portant test for the simulation, and this agreement was 
expected due to the asymptotic character of the Kelbg 
potential as a rigorous weak coupling result. Moreover, 
we observe practically full agreement between MD and 
PIMC results to temperature as low as approximately 
50 000K, corresponding to a coupHng parameter F = 3. 
This is a remarkable extension of "semiclassical" molec- 
ular dynamics into the regime of moderate coupling and 
moderate degeneracy. 

Naturally, below a critical temperature of about 
50 000 K deviations from PIMC results start to grow 
rapidly - the MD results for the energy are becoming 
very low. It is very interesting to analyze the reasons 
for these deviations, as this may suggest directions for 
further improvements. It turns out that the reason for 
these deviations is not a failure of the used quantum pair 
potentials. Thus the only source for the deviations in the 
full simulation can be many-particle effects beyond the 
two-particle level. 

To verify this hypothesis we performed a careful in- 
spection of the microscopic particle configurations in the 
simulation box. At high temperatures, the particle tra- 
jectories are those of a fully ionized classical plasma. At 
temperatures below on Ry, we observe an increasing num- 
ber of electrons undergoing strong defiections on protons 
and eventually performing quasibound trajectories. Most 
of these electrons remain "bound" only for a few classical 
orbits and then leave the proton again. Averaged over 
a long time our simulations are able to reveal the de- 
gree of ionization of the plasma. At the same time we 
observe occasional events of three and more particles be- 
ing at short distances for the duration of one or more 
orbits, refiecting the appearance of hydrogen molecules 
H2, molecular ions H2 etc. 

If the temperature is lowered below approximately 
T = 50 000 K, we observe a strong increase of molecule 
formation and even an aggregation of many molecules 
into clusters with an interparticle distance close to one 
ttB- This turns out to be the reason for the observed 
very low energy because the attractive Coulomb inter- 
action contributions are becoming dominant in the total 
energy. Of course, this behavior is not surprising: while 
all pair interaction processes are modelled correctly even 
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Figure 8: Semiclassical MD results (full lines) for the inter- 
nal energy per hydrogen atom at densities rs= 4, 6 versus 
temperature. The results of restricted PIMC simulations by 
Militzer [4l| are shown for comparison (dashed lines). Sym- 
bols indicate the five temperatures for which MD simulations 
have been performed: T = 31 250, 50 000, 62 500, 125 000 and 
166 670 K (solid lines). The pair approximation breaks down 
around 50 000/S', at the molecule binding energy. 
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Figure 9: Electron-electron and proton-proton radial pair dis- 
tribution functions for a correlated hydrogen plasma with ra = 
4 (left row) and r,= 6 (right row) for T = 125 000, 61 250 and 
31250 K (from top to bottom). 



at low temperature (which is assured by the fit parame- 
ters in the improved Kelbg potentials) , as soon as three or 
more particles are being closely together, three-particle 
and higher order correlations are becoming increasingly 
strong (they, in particular, account for the formation of 
the larger bound state complexes described above) . How- 
ever, it was just the approximation used in the deriva- 
tion of the quantum potentials that three-particle and 
higher correlations can be neglected which was the ba- 
sis for the use of pair potentials in modelling the whole 
A^-particle system. While molecular dynamics, of course, 
includes any level of correlations, the use of the present 
potentials means that quantum corrections to three-body 
(and higher order) interactions are not adequately cap- 
tured. Therefore, it is no surprise that this approxima- 
tion breaks down at sufficiently low temperature, and this 
break down occurs around the temperature correspond- 
ing to the binding energy of hydrogen molecules. From 
this we can conclude that molecule formation sets the 
limit of the applicability of the present semiclassical MD 
simulations. 

Let us now turn to a more detailed analysis of the spa- 
tial configuration of the particles in the MD simulations. 

In Fig. El the radial pair distributions between all par- 
ticle species with the same charge are plotted at two den- 
sities. Consider first the case of T = 125 000 K (upper 
panel). For both densities, all functions look qualitatively 
the same, showing a depletion at zero distance due to 
Coulomb repulsion. Besides, there are differences which 
arise from the spin properties. Electrons with same spin 
show a slightly broader "Coulomb hole" around r = 
than the protons, because the Pauli principle yields an 
additional repulsion of the electrons (this effect is much 



weaker for two protons due to their much smaller de 
Broigle wavelength) . This trend is reversed at lower tem- 
perature, see middle panel, which is due to the formation 
of hydrogen atoms, see also Fier. 1111 below. In this case, 
electrons (their trajectories) are "spread out" around the 
protons giving rise to an increased probability of close 
encounters of two electrons in different atoms compared 
to two protons. 

Now, let us compare electrons with parallel vs. elec- 
trons with antiparallel spins. In all cases, we observe a 
significantly increased probability to find two electrons 
with opposite spin and small distances below one Bohr 
radius which is due to the missing Pauli repulsion in this 
case. This trend increases with lowering of temperature 
due to increasing quantum effects and thus convincingly 
confirms that spin effects are correctly reproduced in our 
MD simulations. 

Before analyzing the lowest temperature in Fig. El let 
us consider the electron-proton pair distributions which 
are shown in Fig. ^1 With lowering of temperature, 
we observe a strong increase of the probability to find an 
electron close to a proton. In contrast to the classical case 
of a collapse (see above), here this probability is finite. 
Multiplying these functions by gives essentially the 
radial probability which is plotted in Fig. El Here, low- 
ering of temperature leads towards formation of shoulder 
around l.Sas which is due to the formation of hydro- 
gen atoms. This conclusion is confirmed by inspection 
of the corresponding quasibound electron trajectories as 
discussed above. At this temperature, the observed most 
probable electron distance is not las as in the hydrogen 
ground state which is due to the considerable kinetic en- 
ergy of the particles leading to a larger average radius of 
the classical quasiclosed orbits We expect that at lower 
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Figure 10: Electron-proton radial pair distribution functions 
at rs= 4 (left figure) and rs= 6 (right figure) and five tem- 
peratures: T = 166 667, 125 000, 62 500 and 50 000 K. 



temperature, the most probable radius would tend to- 
wards laB, but this temperature range is not realistically 
modelled due to molecule and cluster formation. 




Figure 11: Electron-proton radial pair distribution functions 
multiplied by r . Same data as in Fig.lTnl 

While the description of correlated complexes of more 
than two particles is certainly beyond the present pair 
approximation model, nevertheless, several features of 
a partially ionized and partially dissociated hydrogen 
plasma are reproduced correctly. At 62 500 K and r^. = 6 
(right center part of Fig. |2l the simulations show a first 
weak signature of molecule formation - see the maximum 
of the p-p pair distribution function around r = 2aB and 
the maximum of the pair distribution function of elec- 
trons with antiparallel spins around r = l.Sos. Further 
lowering of the temperature by a factor of two (lower 
panel of Fig. |2l confirms this trend: the p-p functions 
exhibit a clear peak very close to r = IAub - the theo- 



retical p-p separation in iJ2-molecules. At the same time, 
also the e-e functions have a clear peak around r = O.bas, 
in the case of opposite spins, and r — 1.2aB, for parallel 
spin projections. The first case comes rather close to the 
true quantum mechanical H — H bound state (singlet) 
with the electron wave function predominantly concen- 
trated between the two protons. On the other hand, this 
electron peak should also extend to the right of the p-p 
peak, and no such pronounced peak would be expected 
for electrons with the same spin. 

Nevertheless, we may conclude that even formation 
and spatial dimension of hydrogen molecules appear to be 
captured surprisingly well in these simulations. The main 
difficulty appears to arise not on the level of four-particle 
correlations but on the level of six particle correlations: 
in the simulations nothing prevents two "bound" atoms 
from binding to a third and more atoms. The overall at- 
tractive Coulomb interaction makes it, below 50 000 K, 
energetically favorable to form large clusters consisting 
of more than two atoms, explaining the strong decrease 
of the internal energy at the T = 31 250 K, cf. Fig. |H| 
In reality, complexes of two molecules do exist, but they 
have a very low binding energy which is due to a subtle 
compensation effects arising from repulsive exchange in- 
teraction between the electrons which go far beyond the 
level of pair interactions . 



VI. USE OF THE QUANTUM PAIR 
POTENTIALS IN DENSITY FUNCTIONAL 
THEORY 

The effective quantum potentials have been introduced 
to represent the equilibrium two particle density matrix 
and subsequently generalized to incorporate many-body 
Coulomb coupHng effects. There are other many body 
coupling effects due to degeneracy or exchange correla- 
tions. For some applications, it may be useful to incor- 
porate these directly in the effective pair potentials to 
extend their validity to still lower temperatures, as was 
demonstrated on the example of classical MD above. 

In this section, we describe the usefulness of the effec- 
tive quantum potentials for a completely different theo- 
retical approach - density functional theory (DFT). In 
doing so, the role of effective quantum potentials with 
degeneracy effects is illustrated as well. 

DFT is a formal structure in which non-perturbative 
approximations can be introduced to describe strong 
coupHng effects Although there are both classical 

and quantum versions of DFT, the classical form does 
not apply to a system of electrons and positive ions due 
to the Coulomb divergence. One possibility is to postu- 
late a classical statistical mechanics using the effective 
quantum potentials described above which allows to re- 
move the singularity. Alternatively, the proper quantum 
formulation can be used from the outset and the effective 
quantum potentials "derived" as a tool in the process of 
computing properties of interest 01. This second ap- 
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proach will be used here. 

In essence, DFT is a variational means to derive an 
equation for the charge density induced by an external 
potential. If that potential is taken to be the same as 
the potential of one of the charges in the system, the 
resulting density is in fact formally identical to the equi- 
librium pair correlation function, or diagonal element of 
the two particle density matrix. The density obeys a 
known nonlinear integral equation - a generalization of 
the Boltzmann-Poisson equation. However, in practice, 
the direct solution this equation is seldom attempted. 
Instead an equivalent set of self-consistent one particle 
Schrodinger equations, the Kohn-Sham equations [45| . 
are solved first to construct the charge density. Yet it 
might be very useful to recall the existence of an alter- 
native direct approach which becomes practical if an ap- 
propriate quantum pair potential is introduced. This is 
illustrated in more detail as follows. 

Consider a quantum system in the presence of external 
sources that can be described by an additive potential 



V 



(41) 



1=1 



Here a denotes a species and q.io, is the position operator 
of particle i of species a. The caret on the potentials is 
used to distinguish the quantum operator from its cor- 
responding function. In general, each species may have 
a different form for the coupling to the external sources. 
The potential also can be written in terms of the density 
operators for each species 



I dr}_^y„(r)n„(r), 



n„(r) 



No, 

E 

i=l 



(42) 

The details of the remainder of the Hamiltonian are not 
important at this point. For this many-body system with 
external sources the theorems of density functional the- 
ory apply in the following form. First, a functional of the 
average densities, ^^(r), averaged over an equilibrium 
grand canonical ensemble is constructed (the generaliza- 
tion to other equilibrium ensembles has been carried out). 
This is done in two steps. First, the equilibrium grand 
potential for the system is considered formally 



PVLe = - In ^ Tr ( 



-/3(h-i:. 







(43) 



The density for the various species is obtained (formally) 
by functional differentiation with respect to the poten- 
tials 



fie = Vie ({^a " ^q}) , "^ea {r) 



(44) 



The density equation is inverted (formally) to get the 
external potentials as functionals of the average densities 



and a Legendre transformation is performed to construct 
the free energy as a functional of the densities rather than 
the chemical potentials 

F({nea}) - n,{{^i^-Vo?s) (46) 

+ dr[^a- Va (r I {necr})] 

a '' 

The crucial second step is to extend this functional to 
arbitrary density fields 



F{M) Fi{n}). 



(47) 



The main task of density functional theory is now to 
construct the density functional 

nv {{n}) = F{{n}) - Jdr (^i, - (r)) n„ (r) , (48) 

where, in this definition, Va (r) is not considered to be 
a functional of the density. The main theorem of den- 
sity functional theory is then that this functional has an 
extremum at the equilibrium density 



Sny {{n}) 
Sn 



= = 



SF{{n}) 
Sn 



- [^J.a - Va (r)] 



(49) 



Va=Vair I {rw}), 



(45) 



Furthermore the value of the functional at the equilib- 
rium density is clearly the equilibrium grand potential 

nv{{n}) = n{{t^a-Va}). 

In practice, an approximate free energy functional 
F {{n}) is written and Eq. (I49II is solved to obtain the 
equilibrium density. This density is then used to eval- 
uate the equilibrium grand potential and determine all 
equilibrium thermodynamic properties. Structural prop- 
erties can be obtained as well by choosing the external 
potential at the end to be the same as that for interaction 
among the system particles. In other words, the source 
is chosen to be a particle of the same type as those com- 
prising the many-body system. The densities Uea become 
equilibrium pair correlation functions. 

How should the functional F{{n}) be constructed? 
There is clearly a part associated with an ideal gas, and 
an energy due to the direct Coulomb interactions. These 
can be identified explicitly. In addition there are the more 
difficult parts due to exchange and correlations. Conse- 
quently, it has become standard practice to write the free 
energy as 

F[n] = ^^(0)({n})+ (50) 
drdr'VaAr - r')nair)n,(r') + F,,{{n}), 

a, a 

where F^^'' ({"■}) is the free energy for the non-interacting 
system, the second term is the contribution from the di- 
rect Coulomb interaction, and Fj.c{{n}) denotes the re- 
maining contributions due to interactions from exchange 
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and correlations. Then the extremum condition H4t)|l be- 
comes 01 

Vi°Hr\M) = V^ir)+ (51) 



J2 j dvdv'V^,{v-v')n„ (r') 



5na (r) 



with Va^ (r I n) denoting the functional Ij45|l for the ideal 
gas. Determination of this functional is the central issue 
of the discussion here, and we will show that it is closely 
related to the Kelbg potential analyzed in the bulk of this 
paper. 

The definition of the functional vi°^ (r | n) is straight- 
forward from the representation of the density for an ideal 
Fermi gas in the external potentials 



na{r) = {r\ e' 



-1 

1 \ \r) 



(52) 



This is a single particle problem. The right side is clearly 
a functional of Va through the dependence of the eigen- 
values of — I- Va on the form of the external poten- 
tial. Interestingly, even at the level of the ideal gas de- 
termination of this functional is non-trivial. In the non- 
degenerate limit this equation for the density becomes 



Ua (r) 



r e 



(53) 



If the external potential is chosen to be a Coulomb source, 
then Ij53|l becomes equivalent to the diagonal elements of 
the two particle density matrix in relative coordinates 
which has exactly the form of the pair distribution func- 
tion used to define the effective quantum pair potential, 
cf. Eq. 

Once the exchange and correlation free energy func- 
tional is specified (guessed), llfiTll provides a set of closed 
classical integral equations for the equilibrium densities. 
As will be seen below, a leading approximation is the 
usual Boltzmann-Poisson representation in terms of semi- 
classical potentials. The primary technical difficulty in 
this prescription is the determination of Va\T^ \ n). 
Kohn and Sham noted that Il5ip defines an effective sin- 
gle particle potential and therefore is formally equivalent 
to the ideal gas in this effective potential. Therefore, the 
solution can be constructed by solving the one particle 
Schrodinger equation whose potential is the right side of 
(EH, and calculating the densities from the associated 
form H52|l self-consistently 



fia {r) 



E 



\^,{r)\\ (54) 



This avoids the difficult problem of finding the functional 
Va'' (r I n) but at the cost of having to solve a set of self- 
consistent Schrodinger equations. 



Consider instead an approximate evaluation of the po- 
tential vi°^ (r I n) in terms of an effective quantum po- 
tential Ua{r) defined by 



(55) 



The first equality is similar to a finite temperature 
Thomas-Fermi representation, with a local chemical po- 
tential given by /ia (r) = fia — Ua{r). An important 
difference discussed below is that Ua{r) ^ Va [r). The 
functional relationship of ria (r) to ^a {t) and hence to 
Ua{r) is that for an ideal gas and is well-known. The sec- 
ond equality of Il55|l defines the semi-classical potential 



Ua{r I Vi°^) as a functional of Vq . This relationship of 



(0) 



Ua{r I Va'') to Vq"'' is more difficult to unfold. However, 
it is straightforward to discover it for weak coupling of 
the system to the perturbing potential. The analysis is 
similar to the derivation of the Kelbg potential and will 
not be repeated here. Formally make the replacement 
Va'' — > AVi^"* in II55|I with the corresponding dependence 
on A inherited by Ua{r). Then perform the expansion of 
Un(r) to first order in A, setting A = 1 at the end, to get 



;(0) 



Ua{r) 



dv'TTa {V~V')V^°HV') 



(56) 



where tTq (r, r') is the well-known static Hnear polariza- 
tion function in random phase approximation, 



7r„(r) = (27r)-^ / dr e**^ "" S=« (fc) 



(57) 



~ n\ f Fa{p - hk) - Fa{p) 



(59) 



dua J (27rft)3 p2 _ (p _ ;jk)2 
containing the Fermi distribution 



Fa (p) = e 



1 



In this approximation, the functional relationship be- 
tween the density and the potential is now known 



n.i (r) 



rfp 



J^+J dr'7r(r-r')vf>(r')-pi) 



(60) 

Now it is straightforward to improve this results by sub- 
stitution of H51|l into the right side of IjftOII which gives a 
generalization of the Thomas- Fermi approximation to in- 
clude strong coupHng effects. However, even if Fxc ({"}) 
is neglected the result is the Thomas-Fermi approxima- 
tion in terms of the potential 



^^(r) = / dr'TTa{r-r')Va{r') 



(61) 
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rather than the bare potential Va(r), which has short 
ranged divergences for opposite charge interactions. The 
result here in terms of the nonlocal effective quantum 
potential appears to be a new one that cures some of 
the well-known problems of the "local approximation" 
Thomas-Fermi theory. As indicated below, Va{j^) be- 
comes just the Kelbg potential in the non-degenerate 
limit. The result II6()|I with Ij51|l is a non- linear integral 
equation for the density, including both strong coupling 
and degeneracy effects. There is no longer any need to 
solve the Kohn-Sham equations and the problem is one 
of purely classical analysis. 

It is instructive to consider the non-degenerate Hmit. 
In that case the polarization function is evaluated using 

^0 (p) ~* e ^{^"^o, f^") _ Furthermore, Eq. Il6()ll simplifies 
to 



-/3(7„(r) 



Ua (r) naC 



VW(^') = / dr'7r-i(r-r')C/a(r'). 



(62) 
(63) 



Use of these in the DFT equation 115 111 gives the closed 
equation for the densities 



In 



Ha (r) 



-l3Vo. (r) - j dvdv' Va.a{Y - r')n„ (r') 



The potentials Va (r) and Vaa{Y — r') are "regularized" 
by the polarization function, e.g.. 



(r) = / dr'^„ (r - r') 14 (r') . 



(65) 



It is possible to show |4fl| that in this non- degenerate 
limit Va (r) is just the original Kelbg potential, Eq. 1^. 
Therefore, in the weak coupHng limit where Fxc {{n}) can 
be neglected l(6ijl becomes the usual Boltzmann-Poisson 
equation with effective quantum potentials given by the 
Kelbg potential Q. 

In summary, DFT applications can be performed in a 
semi-classical form without solving the Kohn-Sham equa- 
tions by introducing effective quantum potentials. This 
can be done in a weak coupling approximation similar to 
that first described by Kelbg and yields a closed analyti- 
cal result Q. Based on the results of the above analysis, 
it can be expected that this approach can be extended 
by incorporating as well effects of degeneracy by using 
for the density Eq. ll52ll instead of l(62|l . Furthermore by 
using improved quantum pair potentials - along the lines 
of the improved Kelbg potentials discussed in the previ- 
ous sections - an accurate treatment of the pair problem 
is achieved laying the foundation for advancing DFT to 
the regime of strong coupling. 



VII. CONCLUSION 

In this work we presented an analysis of generalized 
quantum pair potentials. Extending the work of Kelbg 
and others we investigated in detail effective off-diagonal 
and diagonal quantum pair potentials for a correlated hy- 
drogen plasma including spin effects. We studied the 
accuracy of these potentials by an extensive comparison 
with the exact solutions of the Bloch equation. Further, 
we proceeded to an analysis of improved diagonal quan- 
tum pair potentials by correcting the value of the Kelbg 
potential at zero particle separation. Excellent agree- 
ment with the exact solutions of the two-particle Bloch 
equations could be achieved with the help of a single 
temperature-dependent fit parameter for which an accu- 
rate analytical Fade formula was presented. This lead 
to significantly improved diagonal pair potentials com- 
pared to the original Kelbg potential. Moreover, these 
potentials are explicitly spin-dependent and retain the 
advantage of a closed analytical expression. 

These potentials have been appHed in path integral 
Monte Carlo and "semiclassical" molecular dynamics sim- 
ulations of dense hydrogen and were found to give accu- 
rate results over a wide range of parameters. One im- 
portant conclusion, of relevance to PIMC simulations, 
is that the off-diagonal potential gives essentially more 
accurate results (or more rapid convergence) than its di- 
agonal limit, quantitative estimates have been provided. 

Furthermore, we have demonstrated that the spin- 
dependent improved diagonal potentials are of high use 
for "semiclassical" molecular dynamics simulations of par- 
tially ionized plasmas. Our analysis revealed that with 
these potentials one can successfully simulate dense hy- 
drogen up to moderate coupling and degeneracy, from the 
fully ionized to the partially ionized regime. The present 
potentials allow us to correctly model dense plasmas up 
to temperatures as low as the molecular binding energy. 
Further improvements are possible, including the descrip- 
tion of molecular hydrogen, but they require to include 
three-particle and four-particle correlations and exchange 
effects beyond the two-particle level. 

Finally an intimate relation of the quantum potentials 
to density functional theory has been explored which al- 
lows for DFT calculations without the need to solve the 
Kohn-Sham equations. 
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